
|
Dr Jon Shiach Senior Lecturer Department of Computing and Mathematics Manchester Metropolitan University Email: j.shiach@mmu.ac.uk |
|
Dr Stephen Lynch Reader Department of Computing and Mathematics Manchester Metropolitan University Email: s.lynch@mmu.ac.uk |
|
Dr Killian O'Brien Senior Lecturer Department of Computing and Mathematics Manchester Metropolitan University Email: k.m.obrien@mmu.ac.uk |
|
Click here to download the files
Chaos is a mathematical concept where small changes a system can result in very large differences in the behaviour of the system. This was noticed by Edward Lorenz in 1961 when he was using mathematics to model weather. He wanted to repeat a simulation, but since the solutions took a long time for his computer to calculate Lorenz started the second run from halfway through using solutions from the first run. Instead of the second simulation matching the first, the solution soon diverged producing wildly different weather predictions. Mathematicians are used to solutions to equations being predictable, the solution to equations should not change no matter how many times you solve them so this unpredictability is concerning. Lorenz realised that this was because he entered the solutions from the first run correct to 3 decimal places as opposed to the 6 decimal places the computer used to store the numbers. This slight change, which was less that 0.1 percent, caused the solutions to be very different. Lorenz presented his discovery in a talk titled "Predictability; Does the flop of a butterfly's wings in Brazil set off a tornado in Texas?" which has led to this phenomenon being known as the butterfly effect.
This Jupyter notebook introduces readers to chaos theory through the mathematical modelling of a double pendulum and the motion of bodies acting under gravitation. Whilst some of the mathematics required is ordinarily taught at degree level, to follow the content of this notebook you just need an understanding of differentiation and a little curiosity. This notebook contains some Python code to perform calculations and produce animations. No programming experience is necessary to follow the examples, however, readers can make changes to the code and execute to see what affects any changes has on the output.
Jupyter notebooks are documents that combine text and Python code which allow readable documents such as this one to be created that contain executable code used to perform calculations. To run code in a notebook simply click on the code and then on the run button, or alternatively, press the ctrl + enter keys at the same time. Since we will be using commands to perform calculations and producing animations we need to import some commands to help us to do this. Run the code below to import the commands.
from animations import *
Now we can perform some calculations. The Python code below calculates the length of the opposite side of the right-angled triangle from the example calculation above and prints the result. Note how the equation $\textsf{opposite} = 5 \sin(0.6435)$ is entered in a similar way to how we write it on a piece of paper. Can you add a couple of lines of code to calculate the length of the adjacent side as well?
opposite = 5 * sin(0.6435)
print(opposite)
2.999995564825018
Before we look at examples of chaotic behaviour we are first going to introduce the concept of differentiation. Consider the graph of the curve of $y = f(x)$ where $f$ is some function of $x$. Lets say we want to calculate the gradient (or steepness) of the curve $y = f(x)$ at the point $x_0$ denoted here by the red line. The gradient of a line is calculated by dividing the change in the vertical direction by the change in the horizontal direction, i.e.,
$$\textsf{gradient} = \frac{\textsf{change in $y$}}{\textsf{change in $x$}}.$$If we represent the change in $y$ by $\Delta y$ and the change in $x$ by $\Delta x$ then the gradient of the green line is
$$\frac{\Delta y}{\Delta x} = \frac{f(x_0 + \Delta x) - f(x_0)}{\Delta x}.$$Of course looking at the diagram we can clearly see that the gradients of red and green lines are different. If we reduce the value of $\Delta x$ so that the two points $x_0$ and $x_0 + \Delta x$ are very close together then the gradients of the two lines will be very similar. We can write this mathematically as
$$\frac{\Delta y}{\Delta x} = \lim_{\Delta x \to 0} \frac{f(x_0 + \Delta x) - f(x_0)}{\Delta x}.$$This means "what is the limit of the fraction as $\Delta x$ gets close to zero?". To simplify the presentation we let $\mathrm{d}x$ and $\mathrm{d}y$ represent infinitely small changes in the values of $x$ and $y$ ($\mathrm{d}$ is used because $\Delta$ is the Greek character delta which corresponds to the Latin character d) so we write
$$\text{gradient} = \frac{\mathrm{d}y}{\mathrm{d}x}.$$$\mathrm{d}x$ and $\mathrm{d}y$ are known as differentials), $\dfrac{\mathrm{d}y}{\mathrm{d}x}$ is known as the derivative of $y$ and the process of determining the derivative is known as differentiation.
Differentials are a way of determining the rate of change of one variable with respect to another. Consider the velocity of a moving object
$$\text{velocity} = \frac{\text{distance travelled}}{\text{time taken}}.$$So we can think of velocity, $\mathbf{v}$, as the change in the position, $\mathbf{r}$, with respect to time, $t$, i.e.,
$$\mathbf{v} = \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}.$$Acceleration, $\mathbf{a}$, is the rate of change of velocity with respect to time
$$\mathbf{a} = \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}.$$Since we know that velocity is the rate of change in the position with respect to time then we can write acceleration as
$$\mathbf{a} = \frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} \right) = \frac{\mathrm{d}^2\mathbf{r}}{\mathrm{d}t^2},$$which is known as the second-order derivative of $\mathbf{r}$ with respect to $t$.
A differential equation is an equation that is written using a function of the variables and their derivatives. Differential equations allow us to model all sorts of things such as the motion of planetary bodies in the solar system, the concentrations of elements in a chemical reaction and fluctuations in the stock market.
For example, consider Newton's second law of motion which states that "the force acting on a body is equal to its mass times its acceleration", i.e., $F = m\,\mathbf{a}$. We can write this as a differential equation
\begin{align*} F &= m \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} & & \text{or} & \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} &= \frac{F}{m}. \end{align*}So if we know the force acting on a body and the mass of the body we can solve this differential equation to give us the velocity and therefore we know how the body will move.
The problem with solving differential equations is that most of them cannot be solved analytically, in other words, we cannot solve them using a pen and paper. Instead we use computers to calculate a good estimate of the solution. Here we will be using the programming language Python to perform the calculations required to solve differential equations.
We can use derivatives to produce a mathematical model of simple pendulum which is shown in the diagram here. A pendulum consists of weight that is suspended from a fixed point known as a pivot. When the pendulum is displaced to one side of the pivot and released, the force of gravity acting on the weight causes it to move back of forth underneath the pivot.
Here $m$ is the mass of the weight, $L$ is the length of the rod that connects the weight to the pivot, $\theta$ is the angle between the rod and the vertical and $g$ is the acceleration due to gravity (on Earth this is approximately equal to 9.81ms$^{-2}$). Using Newton's second law of motion we can derive the following differential equations that model the motion of the pendulum
\begin{align*} \frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} &= -\frac{g}{L} \sin(\theta), && (\textsf{the acceleration of the angle $\theta$}). \end{align*}The code below defines the differential equations that model the pendulum in the function pendulum as well as the length and displacement angle before solving the differential equations and plotting the solution.
# Define the differential equations that models a pendulum
def pendulum(t, y):
return y[1], -g / L * sin(y[0])
# Define pendulum parmeters
L = 1 # length of the pendulum
theta = 1 # initial displacement angle (in radians)
g = 9.81 # acceleration due to gravity
tmax = 5 # max time for the simulation
fps = 30 # frames per second for the simulation
# Solve the differential equations
t = linspace(0, tmax, tmax * fps + 1)
y0 = [theta, 0]
sol = solve_ivp(pendulum, [0, tmax], y0, t_eval=t)
# Calculate (x,y) co-ordinates of the pendulum
x, y = zeros((2, len(t))), zeros((2, len(t)))
x[1,:], y[1,:] = L * sin(sol.y[0,:]), - L * cos(sol.y[0,:])
# Plot the solution
anim = plotpendulum(t, x, y)
HTML(anim.to_jshtml())
IntProgress(value=0, description='Animating', max=151)
Note that the pendulum oscillates back and forth in a regular motion and returns to the same initial starting point. Since there are no other forces such as friction or air resistance acting on the pendulum it will continue to move in this way forever. This solution is predictable and non-chaotic.
We have seen that a single pendulum behaves in a regular and predictable way but what happens if we connect another pendulum to our pendulum? Consider the diagram of a double pendulum, the first pendulum has a weight of mass $m_1$ a rod length of $L_1$ and initial displacement angle $\theta_1$, the second pendulum, which is connected to the weight of the first pendulum, and has a weight of mass $m_2$ a rod length of $L_2$ and initial displacement angle $\theta_2$.
The differential equations that model the motion of the double pendulum as a bit more complicated than for the single pendulum so are not given here (inquisitive readers may wish to see them here). The code below solves the differential equations that models a double pendulum and produces an animation showing the motion of the pendulum.
# Define the differential equations that models a double pendulum
def doublependulum(t, y):
t1, t2, dt1, dt2 = y
c, s = cos(t1 - t2), sin(t1 - t2)
dt3 = (m[1] * g * sin(t2) * c - m[1] * s * (L[0] * dt1 ** 2 * c + L[1] * dt2 ** 2) \
- g * (m[0] + m[1]) * sin(t1)) / L[0] / (m[0] + m[1] * s ** 2)
dt4 = ((m[0] + m[1]) * (L[0] * dt1 ** 2 * s - g * sin(t2) + g * sin(t1) * c) \
+ m[1] * L[1] * dt2 ** 2 * s * c) / L[1] / (m[0] + m[1] * s ** 2)
return dt1, dt2, dt3, dt4
# Define pendulum and gravity parmeters
N = 1 # number of double pendulums
theta = [2, 2] # initial displacement angles (in radians)
L = [1, 1] # lengths of the pendulums
m = [1, 1] # masses of the pendulums
tmax = 10 # max time for the simulation
fps = 60 # frames per second for the simulation
# Solve differential equation
t = linspace(0, tmax, tmax * fps + 1)
y0 = theta + [0, 0]
sol = solve_ivp(doublependulum, [0, tmax], y0, t_eval=t)
# Calculate the (x,y) co-ordinates of the pendulum
x, y = zeros((3 * N, len(t))), zeros((3 * N, len(t)))
for i in range(2):
x[i+1,:] = x[i,:] + L[i] * sin(sol.y[i,:])
y[i+1,:] = y[i,:] - L[i] * cos(sol.y[i,:])
# Plot solution
anim = plotdoublependulum(t, x, y)
HTML(anim.to_jshtml())
IntProgress(value=0, description='Animating', max=601)
Here we no longer have the regular periodic motion that we saw with the single pendulum, instead the motion of the second weight follows a complicated path. Try experimenting with the displacement angles, masses and lengths to see what affects this has on the motion of the pendulum.
What you may notice about the double pendulum model above is that changing the initial parameter cause big changes to the motion of the pendulum. This can be seen more clearly below where two near idential double pendulums are modelled at the same time. The only difference between the two pendulums is that the second pendulum has the displacement angle $\theta_2$ changed by a very small amount (here the difference is diff = 0.001 which is equivalent to 1 thousandth of a radian).
# Define pendulum and gravity parameters
diff = 0.001 # difference in the θ_2 angle between each run
theta = [2, 2] # initial displacement angles (in radians)
L = [1, 1] # lengths of the pendulums
m = [1, 1] # masses of pendulums
tmax = 10 # max time for the simulation
fps = 60 # frames per second for the simulation
# Solve differential equation
t = linspace(0, tmax, tmax * fps + 1)
y0 = theta + [0, 0]
x, y = zeros((6, len(t))), zeros((6, len(t)))
for i in range(2):
sol = solve_ivp(doublependulum, [0, tmax], y0, t_eval=t)
y0[1] += diff
# Calculate the (x,y) co-ordinates
for j in range(2):
x[3*i+j+1,:] = x[3*i+j,:] + L[j] * sin(sol.y[j,:])
y[3*i+j+1,:] = y[3*i+j,:] - L[j] * cos(sol.y[j,:])
# Plot double pendulum
anim = plotdoublependulum(t, x, y)
HTML(anim.to_jshtml())
IntProgress(value=0, description='Animating', max=601)